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The linear instability of Lotka-Volterra orbits in the homogenous manifold of a two-patch system 
is analyzed. The origin of these orbits instability in the absence of prey migration is revealed 
to be the dependence of the angular velocity on the azimuthal angle; in particular, the system 
desynchronizes at the exit from the slow part of the trajectory. Using this insight, an analogous 
model of a two coupled oscillator is presented and shown to yield the same type of linear instability. 
This unables one to incorporate the linear instability within a recently presented general framework 

\^ • that allows for comparison of all known stabilization mechanisms and for simple classification of 

f*"*) ' observed oscillations. 
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Population oscillations in prey predator systems, as predicted by the Lotka-Volterra equations are known to 
(-^ . be unstable with respect to additive and multiplicative noise. This instability must lead to the extinction of one 
£S1 ' of the interacting species, a fact that has been confirmed in various experiments for well-mixed populations (2j. 
t— I , The persistence of natural prey-predator and host-parasitoid systems, thus, is commonly attributed to their spatial 
structure, such that migration between desynchronized patches yields an inward flow toward the coexistence fixed 
point and is responsible for the sustainability of the oscillations Q. In fact, spatially extended systems tend to support 
finite amplitude oscillation [4|]. The stabilization of such oscillations is considered to be a major factor affecting species 
conservation and ecological balance @. 

The main challenge, thus, is to understand the conditions for the appearance of desynchronization in diffusively 
coupled patches, since diffusion tends to synchronize these patches so after a while the whole system flows to the well 
mixed, unstable limit @. One of the solutions to that problem was presented by Jansen @, S 0- It turns out that 
i the trajectories far from the fixed point become unstable if the inter-patch migration rate of the predator is much 

larger than that of the prey. Jansen used Floquet analysis to show that instability. In this paper, we try to elucidate 
t-H . the underlying mechanism that yields Jansen's instability, to generalize it in the framework of the recently presented 
coupled oscillator model, and to discuss the conditions under which one may observe the stabilizing effect of Jansen's 
mechanism, like oscillation amplitude that grows under noise until it reaches the first unstable orbit. 

The Lotka-Volterra predator-prey system is a paradigmatic model for oscillations in population dynamics It 
, describes the time evolution of two interacting populations: a prey (6) population that grows with a constant birth 
__ • rate a in the absence of a predator (the energy resources consumed by the prey are assumed to be inexhaustible), 
\Q \ while the predator population (a) decays (with death rate /x), in the absence of a prey. Upon encounter, the predator 
may consume the prey with a certain probability. Following a consumption event, the predator population grows and 
the prey population decreases. For a well-mixed population, the corresponding partial differential equations are, 
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— = -fia + Xiab (1) 



at = CT6 ~ Aaa& ' 



where Ai and A2 are the relative increase (decrease) of the predator (prey) populations due to the interaction between 
species, correspondingly. 

The system admits two unstable fixed points: the absorbing state a = b = and the state a — 0, b = 00. There 
is one marginally stable fixed point at a = CT/A2, b = £t/Ai. Local stability analysis yields the eigenvalues ±i^/fia 
for the stability matrix. Moreover, even beyond the linear regime there is neither convergence nor repulsion. Using 
logarithmic variables z = ln(a), q — ln(b) eqs. ((TJ) take the canonical form z = dH/dq, q = —dH/dz, where the 
conserved quantity H (in the ab representation) is, 

H = Aib + A 2 a - [i ln(a) - a ln(b). (2) 

The phase space, thus, is segregated into a collection of nested one-dimensional trajectories, where each one is 
characterized by a different value of H, as illustrated in Figure 3. Given a line connecting the fixed point to one of the 
"walls" (e.g., the dashed line in the phase space portrait, Figure 3), H is a monotonic function on that line, taking 
its minimum H m i n at the marginally stable fixed point (center) and diverging on the wall. Without loss of generality, 
we employ hereon the symmetric parameters /i = er = A 1 = X 2 = 1. The corresponding phase space, along with the 
dependence of H on the distance from the center and a plot of the oscillation period vs. H, are represented in Figure 
3). 



2 




1 2 3 4 0.5 1 

a ' 



FIG. 1: The Lotka-Volterra phase space (left panel) admits a marginally stable fixed point surrounded by close trajectories 
(three of these are plotted). Each trajectory corresponds to single H denned in Eq. p). where H increases monotonically along 
the (dashed) line connecting the center with the a = wall, as shown in the lower right panel. In the upper right panel, the 
period of a cycle T is plotted against H, and is shown to increase almost linearly from its initial value T = 2n / \fji5 close to 
the center. 




o.o-| — , — i — r— ^ — , — , — i I — i — i — , — '] r j 

10000 20000 30000 40000 50000 60000 70000 BO000 

rat 



FIG. 2: The survival probability Q(t) is plotted versus time for a single-patch, noisy LV system. Eqs. |T} (with the symmetric 
parameters) were integrated numerically (Euler integration with time step 0.001), where the initial conditions are at the fixed 
point a = b = 1. At each time step, a small random number rj(t)At was added to each population density, where n(t) £ [—A, A]. 
A typical phase space trajectory, for A = 0.5, is shown in the inset. The system "dies" when the trajectory hits the walls a = 
or b — 0. Using 300 different noise histories, the survival probability is shown here for A = 0.5 (full line), A = 0.3 (dotted 
line) and A = 0.25 (dashed line). Clearly, the survival probability decays exponentially at long times, Q(t) ~ exp(—t/r). As 
expected for a random walk with absorbing boundary conditions, 1/r scales with A 2 . 

Given the integrability of that system, the effect of noise is quite trivial: if a and b randomly fluctuate in time 
(e.g., by adding or subtracting small amounts of population during each time step), the system wanders between 
trajectories, thus performing some sort of random walk in H with "repelling boundary conditions" at H m i n and 
"absorbing boundary conditions" on the walls (as negative densities are meaningless, the "death" of the system is 
declared when the trajectory hits the zero population state for one of the species). This result was emphasized by 
Gillespie [nj for the important case where intrinsic stochastic fluctuations are induced by the discrete character of 
the reactants. In that case, the noise is multiplicative (proportional to the number of particles), and the system flows 
away from the center and eventually hits one of the absorbing states at 0,0 or 0,oo. The corresponding situation for 
a single patch Lotka-Volterra system with additive noise is demonstrated in Figure 4, where the survival probability 
Qit) (the probability that a trajectory does not hit the absorbing walls within time t) is shown for different noise 
amplitudes. 

The Lotka-Volterra system on spatial domains has been investigated, usually in a form of diffusively coupled patches, 
during the last decades. Any patch is assumed to be well mixed, and the flow of the reactants from one patch to its 
neighbors is governed by the density gradient. Clearly, any system of that type, independent of its spatial topology 
(either regular lattice of some dimensionality or some sort of network without isolated nodes) admits an infinite 
number of solutions that correspond to synchronous oscillations of the whole system along one of the H trajectories, 
where the diffusion has no role as there are no population gradients. The simplest example is the two-patch system, 
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described by: 



~dt 

da 2 

~dt 
dbi 

~dt 

dh 
dt 



-fiai + Aiai&i + D a (a 2 - ai) (3) 
— fia 2 + Aia 2 6 2 + D a (a\ — a 2 ) 
(T&i - \ 2 aibi + D b (b 2 - b x ) 
ab 2 - A 2 a 2 6 2 + D b (bi - b 2 ). 



Here the invariant manifold is the two dimensional subspace a\ — a 2l b\ — b 2 . The diffusion, of course, suppresses 
fluctuations and stabilizes the invariant manifold; one may expect, thus, that the single-patch dynamics also capture 
the main features of the extended system, and that the system behaves like a random walker in the invariant manifold 
(with a rescaled noise) and hits the absorbing walls after some characteristic time r, where r scales linearly with the 
noise strength A 2 . 

As a first hint for a stabilizing mechanism, let us consider the total H , 

H T = Hi + H 2 = ai + a 2 + bi + b 2 - ln(aia 2 bib 2 ) . (4) 
With the deterministic dynamics ([3]), Ht is a monotonously decreasing quantity in the non- negative population regime: 



dt \ a\a 2 J \ bib; 

Accordingly, if an orbit on the invariant manifold becomes unstable, the flow will be inward and the population 
oscillations stabilizes. 

While if D a = D b the stability properties of an orbit on the invariant manifold are identical to the stability 
properties of the corresponding single-patch orbit [11] , if the diffusion of both species is different, there is a possibility 
for unstable orbits on the homogenous plane. This option was materialized by Jansen [?}, who considered the set of 
Equations [3] in the limit Db = 0, so that only the predator undergos diffusion. With the transformation: 

ai + a 2 „ 6i + &9 
A = B = (6) 

2 2 W 

ai - a 2 bi-b 2 

o = a = . 

2 2 ' 

one recognizes the homogenous AB manifold and that the 5 9 coordinates measure the deviation from that manifold 
(the heterogeneity of the population). In these coordinates the system satisfies, 

8A 

^- = -fxA + XiAB + XiSO (7) 
ot 

— = aB - \ 2 AB + \ 2 S9 
dt 

_ = -nd + X X A6 + XiB5 - 2D a S 
ot 

f)9 

— = aB - X 2 A9- X 2 B5-2D b 9. 
ot 

Linearizing around the homogenous manifold, The AB dynamic is equivalent to that of a single patch, 

A = -fxA + XiAB (8) 
B = oB- X 2 AB 



and the 6 — 9 linearized dynamic is 

d_fS\_f- ( i + X 1 B-2D a XtA . . .. . 

dt { ) \ -X 2 B a- X 2 A-2D b M n J • 
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FIG. 3: Stability diagram for phase space orbits (ordered by their conserved quantity H) for different values of predator 
diffusion D a , where Db=0. 
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FIG. 4: The angular velocity along some orbits of the Lotka-Volterra dynamic. Fast regions marked in red, slow regions are 
blue. Clearly, the dynamics is slowest when the populations of both species are diluted, and fastest along the dense region in 
the upper-right "shoulder." Note that the velocity gradient along an orbit increases with H. 



One may thus calculate the eigenvalues of the Floquet operator for one period along an orbit of the homogenous 
manifold ([5]). The resulting stability diagram, first obtained by 0], is shown in Figure [3) 

Our first mission is to intuitively explain Jansen's results. First, we notice that the angular velocity along a single 
Lotka-Volterra orbit is not fixed. Figure 0] emphasizes the angular velocity gradient along an orbit. While the inner 
trajectories (close to the fixed point) are almost harmonic with constant angular velocity, the eccentric large H orbits 
admit large variations. In particular, the motion in the dilute population region [close to the unstable empty fixed 
point (0, 0)] is very slow, while in the dense population region the angular velocity is large. 

Following the caricature of an orbit in Fig. we can explain the source of the instability. For a two-patch system, 
if one patch is at point A along the orbit and the other patch at B, since the A patch is moving faster along the line 
it will get closer and closer to B during their flow toward the slow region. The diffusion of the prey plays no role 
along this branch, since the prey density is almost equal, while the predator diffusion may only strengthen that effect. 
Thus, the two patches must (almost) synchronize along this branch. 

The situation is completely different in the exit from the slow region. The patch at D moves much faster than that 
at C, so they will desynchronize. As the predator density along this branch is almost constant, the only factor that 
may avoid desynchronization is the prey migration. In the absence of prey migration, the two patches reach the points 
C and D', where the predator migration produces an inward flow. Figure [HI is now well understood: the inward flow 
happens when the desynchronization interferes with the predator diffusion, as explained. 

Let us consider, now, Jansen's instability in the framework of the coupled nonlinear oscillator toy model, recently 
presented fl2| as a generic tool for the investigation of oscillation stability in diffusively coupled metapopulations. 
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FIG. 5: An orbit of the LV dynamics and its fast and slow regions. As explained in the text, with no prey migration the 
two patches desynchronize in the CD region, thus predator diffusion causes a flow toward the fixed point and stabilizes the 
oscillations. 
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FIG. 6: Phase portrait of the inward flow in the homogenous manifold (average prey density vs. average predator density) for 
two-patch LV system with no prey diffusion and D pre dators = 1- Clearly, the inward flow happens in the C" — D' region of 
Figure O where the desynchronization along the CD branch interferes with the predator diffusion. There is almost no inward 
motion along the rest of the orbit. 



With the intuition gathered from the above example, we want to consider diffusively coupled orbits where the angular 

velocity depends on the radial angle and the diffusing species density is changing along the slowing branch. The 
following equations, 

—± = co(e 1 )y 1 +D x (x 2 -x 1 ) (10) 
ot 

dx 2 

— — = w(6 2 )y 2 + D x (xx - x 2 ) 
ot 

^ = -uj(9 1 )x 1 + D y (y 2 - yi) 

^ = -uj(d 2 )x 2 + D y (y 1 -y 2 ), 



will satisfy these conditions for D x — D, D y — and 



7T 

CJ + Wl cos(# - -). (11) 
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Using (i G 1,2) 

rf = xj + yf 6, = arcban(—) (12) 

(xx + yy) ■ (xy - yx) 

r= 9= , 

and 

T = r 2 — r\ R = r-2 + r\ (13) 

cj) = 6 2 -8 1 ^ = e 2 + 9i, 

one finds that the flow in the invariant manifold satisfies: 

R = (14) 



while the linearized equations for the desynchronization amplitude r and the desynchronization angle <fi satisfy: 



d ( <j> \ _( 2w'($/2) - 2D X cos 2 ($/2) 



2/'^f>/0^ n x sin<[> 



R 



dt\r ) \ -DvRsinW -2^sin 2 ($/2) 



(15) 



Using the Floquet operator technique to analyze the stability of an orbit by integrating (|15|) along a close trajectory 
of fT4")) one finds the stability map presented in Figure [5J where the parameter H of Figure [3J is now replaced by w\ , 
which measures the "eccentricity" of the angular velocity along a circular path. Here, two unstable regions appear, 
for large and small D x . 

It is interesting to point out that in the high D x region, the unstable eigenvalue is positive, while in the small 
D x unstable region it takes negative values. The reason for that is the effect of predator migration. If the effect of 
migration is large, in comparison with the intra-patch dynamics, the two patches should admit (almost) the same 
predator density, and the corresponding points in the 2d phase portrait should stay on the same vertical line (same 
"x" coordinate). The trajectories of the two points representing the patches along an orbit are thus similar to the 
transport of a vertical rod along a circle, keeping the center of the rod on the circular trajectory: the two ends switch 
their role twice (the inner becomes the outer and vice versa) so the rod returns to its original state after a full cycle 
(See Figure [3 left panel). In the low diffusion range, on the other hand, the points remain on a vertical line only 
close to the slow portion of the orbit, where they switch once, but in the fast region they support different predator 
populations, so the "rod" completes its cycle in opposite "phase" (see Figure [3 right panel). 

We now turn to our last point, a comparison of this stabilizing mechanism with the nonlinear, noise induced 
mechanism recently presented by us (l2j . The stability mechanism of (l2| involved with the amplitude dependence 
of the angular velocity (see the upper right panel of Figure [1]) , works as well for system of equal prey and predator 
migration rates and is not based on a linear instability of an orbit. One may ask, thus, how to make a distinction 
between these two mechanisms in real systems. 

In order to make a distinction between amplitude-induced stability [l2j and angular-induced stability Q one should 
compare the corresponding radius of oscillations, where the dominant mechanism corresponds to the smaller radius. 
The amplitude synchronization prediction is that the oscillation radius scales like D/(oj') 2 , where u>' = doj/dr is the 
frequency gradient along the oscillations amplitude (See Figure [TJ upper right panel). This result should be compared 
with the instability radius of [7|, and for small migration rates {D ~ 0.01) it is smaller in few orders of magnitude. 
It seems, thus, that the angular induced instability will be relevant only for relatively large diffusivities, where the 
effect of amplitude-induced instability is suppressed by patches synchronization. 

In order to observe the phase instability, we have simulated the two-patch LV system, where the effect of demographic 
stochasticity, an intrinsic noise that should appear in any system independent of environmental factors, was introduced 
via a noise term proportional to the square root of the population size. The results are presented in Figure [51 for the 
two opposite cases: no predator migration (where one should expect angular instability), and no prey diffusion, where 
no such instability is present. Both cases were simulated for /i = a — 1, so the population (number of individuals in 
each of the species) at the coexistence fixed point is 1/A. The diffusion D — 1 corresponds to the smallest amplitude 
of the last stable orbit and is way too large to allow amplitude desynchronization (see [l2T]). 

Fig. [S] clearly shows the stabilizing effect of angular-induced desynchronization. For small populations at the fixed 
point (large A) the prey diffusion systems reach the absorbing state, while only predator diffusion stabilizes the inner 
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FIG. 7: An illustration of the trajectories of two diffusively-coupled patches, with slightly different initial conditions, projected 
on the invariant manifold. In the strong coupling case, (left panel) the strong predator diffusion forces the two points to be on 
the same vertical line (same predator concentration) along the orbit, hence the phase of the Floquet eigenvalue inverted twice 
along the trajectory, yielding a positive eigenvalue. In the small diffusion limit, the patches posses equal predator density only 
in the slow portion of the orbit, when the intra-patch dynamic is slow with regard to the migration. This leads to trajectories 
like those illustrated in the right panel (points connected by "rod" stand for the population density in equal times), where only 
one sign change happens and the Floquet eigenvalue is positive. 



orbits. Large populations, though, may be stable in both regimes, but the instability cuts the tail of the distribution, 
leaving only a peak close to the "reflecting boundary." 

To conclude, it has been shown that systems where only the predator admits the ability to migrate [a canonical 
examples include herbivore - plant or parasite insect - plant systems, like in the famous example of biological control of 
the Prickly Pear cactus by the moth Cactoblastis cactorum in eastern Australia [HI]] may support sustained oscillations 
in noisy environments. This phenomenon has been explained here, and its cause was traced to the dependence of the 
angular velocity on the azimuthal angle along an orbit lying in the homogenous manifold. This insight allows us to 
incorporate that phenomenon into a generic framework of coupled nonlinear oscillators and to compare that mechanism 
with other stabilizing effects. In a separate publication [lil ]. we intend to put forward a general classification scheme 
for stable population oscillations, a scheme that may be used to typify the observed desynchronization-induced stable 
manifold according to its underlying mechanism. 
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FIG. 8: Stability diagram in the u>i — D x plane for the Floquet operator (same as Figure [3} for the coupled oscillator system 
described by Eqs. (flOjl with D y — 0. The two unstable regions correspond to different signs of the Floquet unstable eigenvalue, 
as explained in the text. 
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FIG. 9: Histograms of the probability density as a function of H, for a two-patch LV system with only prey diffusion (Db = 
1, D a = 0)(b) and only predator diffusion (D a = 1, Db = 0)(a). Both systems were subject to demographic stochasticity, 
modeled by a multiplicative noise proportional to the square root of the population density. In both cases, the probability 
density is concentrated around H = 0; however, Jansen's instability manifests itself in the peak at the instability radius at the 
left panel, caused by the "reflection" from the unstable manifold (note the arrow that indicates the first unstable orbit). The 
log- log plots of that histogram (insets) show that the tail of the distribution is continuous at the right panel, but the probability 
to find the system with H above the instability limit is practically zero. The LV parameters are fx = a = 1 and A = 10 -5 . 
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